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(-H I and semi-smooth Newton methods for the numerical solution of pde constrained optimization with 

control constraints [3l [11] special emphasis has to be taken on the implementation, convergence and 
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^ ■ 1 Introduction and mathematical setting 



X 



We are interested in the numerical treatment of the following control problem 

min(y,«)gyxc/„, J{y,u) := ^||y - ^|li2(o) + f Ikllci 

s.t. ^ ^ 

(1.1) 

—Ay = Bu in fi, 
y = on . 

Here, Q C {d > 1) denotes an open, bounded sufficiently smooth (polyhedral) domain. 
Given some Hilbert space U and some closed, convex admissible set Uad C U for the controls 
and a linear, continuous control operator B : U ^ H^^{Q), the states lie in Y := Hq{^). Let 
us note that also additional state constraints could be included into our problem setting, as 
done in jT] and [2], and also more general (linear) elliptic or parabolic state equations. How- 
ever, all structural issues discussed in the present work are induced by the control constraints, 
hence to keep the exposition as simple as possible state constraints are not considered here. 
Typical configurations of P are 
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Examples. 

m 

(i) U := M™, Y = H^{n), B -.W^ ^ H-\n), Bu := ^ ujFj, Fj £ R-^n), U^a ■= {v G 



aj < Vj < bj}, a, 6 G M"", a < b. 



(ii) U := L^(ri), Y = Hq{U), B = t : L?'{Q.) H ^(^2), t being the canonical injection, 
:= {v £ L'^{n)-a <v<b},a,be L°°{n), a<b. 

Remark 1.1. One may as well consider elliptic equations with Neumann boundary control, 

-Ay + y = in r2, 
dr^y = u on dQ , 

thus setting U := /.^(r), Y = H\Q), U^a ■= {v G ^^(r); a<v<b},a,be L°°(r), a < b. 

Problem P admits a unique solution (y,n) £Yx f/adj and can equivalently be rewritten as 
the optimization problem 

min J{u) (1.2) 

for the reduced functional J(u) := J{y{u), u) = J{SBu, u) over the set Uad, where S : Y* ^ Y 
denotes the (continuous) solution operator associated with —A. We further know that the 
first order necessary (and here also sufficient) optimality conditions take the form 

{J'{u),v - u)u*,u > for all v G (1.3) 

where J'{u) = au + B* S*{SBu — z) = au + B*p, withp := S* {SBu — z) denoting the adjoint 
variable. The function p in our setting satisfies 

— Ap = y — z in $7, 

^ ^ (1.4) 
p = on di}. 

For the numerical treatment of problem (II. ip it is convenient to rewrite (11.30 for o" > 
arbitrary in form of the following non-smooth operator equation; 

^ = Pu., - ^TVJ(n)) Pu^, (^-l^R-^B*p^ , 

with the Riesz isomorphism R : U U* and the gradient V J(m) = R^^J'{u). 

2 Finite element discretization 

To discretize (P) we concentrate on Finite Element approaches and make the following as- 
sumptions. 

Assumption 2.1. 

O C M'^ denotes a polyhedral domain, = U^^, T,- with admissible quasi-uniform sequences of 
partitions {Tj}"^^ of i.e. with hnt '■= maxj diam Tj and Unt '■= minjjsupdiam K;K <Z Tj} 
there holds c < ^ < C uniformly in nt with positive constants < c < C < co independent 
of nt. We abbrevi^ate Th := {Tj}]ii. 
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For G N we set 



Wh ■■= {v £ C70(0); v\^^ G Pfc(T,) for all I < j < nt} =: {cPi, . . . , cP^g), and 

Yh := {r; G W^, v\^^^ = 0} =: (<^i, . . . , </.„) C F, 

n 

with some < n < ng. The resulting Ansatz for yh then is of the form yh = ^ Ui'Pi- Now we 

i=l 

approximate problem (P) by 



s.t. (2.1) 

a{yh,Vh) = {Bu,Vh)Y*,Y for all v/j G y/i, 



where a{y,v) := J^VyVvdx denotes the bilinear form associated with —A. Problem {¥h) 
admits a unique solution {yh,u) £ x [/^d and, as above, can equivalently be rewritten as 
the optimization problem 

min Jh{u) (2.2) 

tlGf/ad 

for the discrete reduced functional Jh{u) '■= J{yh{u)-,u) = J{ShBu,u) over the set f/ad, 
where '■ Y* ^ Y^ C Y denotes the solution operator associated with the finite element 
discretization of —A. The first order necessary (and here also sufficient) optimality conditions 
take the form 

(Jhi^h), V - Uh)u',u > for all v G C/ad (2.3) 

where J^(f ) = av + B* S'f^{ShBv — z) = au + B*ph, with ph := S^{ShBu — z) denoting the 
adjoint variable. The function ph in our setting satisfies 

a{vh,Ph) = iVh - z,Vh)Y',Y for ah Vh G Y^. (2.4) 
Analogously to (II. 3p . for cr > arbitrary, we have 



Pu^^ {-^h - crVJhiuh)) " = " P[/,, {-^R-^B*pj}j . (2.5) 



Remark 2.2. Problem ()2.ip is still infinite-dimensional in that the control space is not 
discretized. This is reflected through the appearance of the projector Pjj^^ in (|2.5p . The 
numerical challenge now consists in designing numerical solution algorithms for problem p.ip 
which are implementable, and which reflect the infinite-dimensional structure of the discrete 
problem [HIS]. 



Next let us investigate the error — n/i||[/ between the solutions u of (11. 2p and Uh of (|2.2p . 
compare [7]. 

Theorem 2.3. Let u denote the unique solution of (II. 2p . and Uh the unique solution of 
(12:2]). Then there holds 



"Ik - UhWu + ^\\y{u) - Vhf < {B*{p{u) -ph{u)),uh - u)u*,u + ^\\y{u) - yh(w)|li2(Q), 

(2.6) 

where Ph{u) := S]l{SBu — z), yh{u) := ShBu, and y{u) := SBu. 
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Proof. Since ()2.2I) is an optimization problem defined on all of t/ad) the unique solution u 
of (jl.2p is an admissible test function in (j2.3p . Let us emphasize, that this is different for 
approaches, where the control space is discretized explictly. In this case we may only expect 
that Uh is an admissible test function for the continuous problem (if ever). So let us test (jl.3p 
with Uh, and (12. 3p with n, and then add the resulting variational inequalities. This leads to 

{a{u - Uh) + B*S*{SBu - z) - B*Sl{ShBuh - z),Uh - u)jj^^^ > 0. 

This inequality is equivalent to 

«lk - ^hWu < {B*{p{u) - ph{u)) + B*{ph{u) - Ph{uh)), Uh - u))u*,u- 

Let us investigate the second addend on the right hand side of this inequality. By definition 
of the adjoint variables there holds 

{B*{ph{u) - Ph{uh), Uh - u)^,^jj = {ph{u) - Ph{uh), B{uh - u))y,y* = 

= aivh - yhiu),ph{u) - Ph{uh)) = j (vhiuh) - yh{u)){y{u) - yh{uh))dx = 

n 

/I I 
{y - yh){y - yh{u))dx < --\\yh - y\\l2(^n) + 2 ~ y/^MIlL^cn) 

SO that the claim of the theorem follows. I 

What can we learn from Theorem 12.61 .'' It tells us that an error estimate for ||u — it/i||{/ is at 
hand, if 

• an error estimate for \\R~^B*{p{u) — Ph{u)\\u is available, and 

• an error estimate for \\y{u) — y/i('u)||L2(Q) is available. 

Remark 2.4. The error ||u — ttft,||{/ between the solution u of problem (II. 2p and Uh of (12. 2p 
is completely determined by the approximation properties of the discrete solution operators 
Sh and St. 



3 Semi-smooth Newton algorithm 

In the following we restrict our considerations to the practically relevant case of the second 
example given in Section [H i.e. we set U = Y = Hq{^1), Uad = {v G L^(r2); a < v < b} 

with a,b £ L°°(il), b — a > a > and o" G M. Also the control operator is the injection 
I : L^{Q) Y*, hence the adjoint B* = i* is the injection from Y into L^{Q,). Below, the 
operators B, B* and R are omitted for notational convenience. The variationally discretized 
problem associated to (P) then reads 



™™{yh,«)eyxL2{n) Jiy,u) := ^\\y - 2||i2(Q) + f Ik 
s.t. 

"h) { a{yh,Vh) = {u,Vh)L2{n) for all -Uh G 
and 

a < u < b, a.e. in . 
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To apply the semi-smooth Newton algorithm proposed in the following, the bounds are re- 
quired to be elements of the finite element space Y^. Let therefore ah,bh G be obtained 
from a, b by interpolation or projection and let us consider the problem 

s.t. 

'hh) S a{yhh,Vh) = {u,Vh)L2{n) for all t;,, G 1",^ 
and 

ah < u < bh, a.e. in Q, 

It is clear that for h > small enough the admissible set < u < b^ is non empty, if 
o-h,bh —5 a,b uniformly, say which can be guaranteed for sufficiently smooth bounds a, 6 
and ah = Iha, bh = Ihb, with Ih denoting the Lagrange interpolation operator or the L^- 
projection. In this case problem (Fhh) admits a unique solution {uhh^Uhh)- Let us assume, 
that this solution exists. 

Lemma 3.1 (Perturbed Bounds). The solutions (uhh^Uhh) and {yh,Uh) of (Phh) o,^d (P/i) 
satisfy the estimate 

\\uhh - Uh\\L2(n) - + ^"'^^"^) ~ o-hh^in) + \\b - bhh^in)) ■ 

Proof. Let 

Then by (j2.5p there holds 

Ikh - '^h\\L^{Q.) < \\ah - a\\L^(n) + Ph - • 
Since is admissible for ¥hh we have 



(3.1) 



{uhh + -S*h{ShUhh - z),ul - Uhh)L^(n) > 



while the definition of gives 



{Uh + -S*h{ShUh - z),Uhh - uDL2(n) > 
since Uhh lies between ah and bh- Adding these inequalities leads to 

'^hh\\l2(^Q^ < ^{SlSh{uh - Uhh)-,Uhh - L'^{n) 

= -{SlShiuh - u\),uhh - Uh)L'2{n) + -{ShSh{ul - Uhh),Uhh - uDi^n) 
and finally we have 

ll^h - '^hh\\l2(^Q) + l^W^hiul - Uhh)\\l2(^n) < ^{ShShiuh - ul),Uhh - uI)l2(Q) 



a 



which combined with (j3.1|) implies the lemma. 
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Now let 

G{v) ■.= v- P[a,b] (^-^P{y{v))^ , and Gh{v) := v - P[a^,b,,] (^-^Ph{yh{v))J , (3.2) 

where for given v G L'^{^1) the functions p,Ph are defined through (jl.4|) and (j2.4p . respectively. 
It follows from the characterization of orthogonal projectors in real Hilbert spaces that the 
unique solutions u,Uh to (jl.ip and (|2.ip are characterized by the equations 

G{u), Gh{uh) = in L\n). (3.3) 

These equations will be shown to be amenable to semi-smooth Newton methods as proposed 
in [3] and [11]. We begin with formulating 

Algorithm 3.2. (Semi-smooth Newton algorithm for (j3.3p ) 

Start with v £ L^(J7) given. Do until convergence 

Choose M G dGh{v). 

Solve M5v = —Gh{v), v := v + Sv. 

If we choose Jacobians M £ dGh{v) with ||M~^ || uniformly bounded throughout the iteration, 
and at the solution Uh the function Gh is SG^-semismooth of order /i, this algorithm is locally 
super convergent of order 1 + Although Algorithm 13.21 works on the infinite dimensional 
space L^(r2), it is possible to implement it numerically, as is shown subsequently. 



3.1 Semismoothness 

To apply the Newton algorithm, we need to confirm that the discretized operator Gh is indeed 
semismooth. To establish this fact we rewrite Gh in the form 

Gh{u) =u - (6- a)P[o,i](^(fe- a)""^(^ - ^{Sl{ShU- z)) - a)j^ + a 

and apply ([llj. Theorem 5.2), with -P[o,i] : M — > M taking the role of ip. Here and in the 
following, for notational convenience we assume a,h G Yh, which is no restriction due to 
Lemma l3. II The smoothing-operator F : ^ L'^ from [11] in our case reads 

F{u) = (6 - a) ( - 1 {SUShU - z)) - a) . 

We note that 

• since we require a,b £ L°°{Q,), b — a > a > with a G M, both (6 — a) and {b — a)~^ are 
in L°°(J7) and the pointwise multiplication by either (b — a) or (6 — a)~^ is a continuous 
endomorphism in L^(r2) for any p. 



the operator F is differentiable with constant derivative for any g > 1. In fact, for 
sufficiently smooth domains the operators Sh and 5^ map L'^{Q,) continuously into 
H'^{Q), which is continuously embedded in L'^{Q) for any q G [1, cxd]. 

P[o^i] : M ^ M is 9i-*[o,i]"Semismooth of order 1, with 



dP[o,i]ix) 



if x ^ [0,1] 

1 ifxG(0,l) . (3.4) 
[0,1] if a; = or a; = 1 
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• for piecewise linear elements the semismooth complementarity condition (5.3) in 
theorem 5.2) holds automatically with 7 = 1. 

Thus we are in the position to apply .theorem 5.2) with a = 1 and qo > r = 2 and 7=1 
and obtain 

Theorem 3.3. The function defined in liS. ^)) is dGh-semismooth of order /i < There 
holds 

dGh{v)w = w + ^^^[a,b] (^-^PhiVhiv))^ ■ {SlShw) , 

where the application of the differential dP[a,b] ^^^^^ multiplication by S^S^w are pointwise 
operations a.e. in fl. 

Remark 3.4. In [5] the mesh independence of the superlinear convergence is stated. Recent 
results from [T2] indicate semismoothness of G of order i as well as mesh independent q- 
superlinear convergence of the Newton algorithm of order |, if for example the modulus of 
the slope of —^p{y{u)) is bounded away from zero on the border of the active set, and if 
the mesh parameter h is reduced appropriately. This is the key to our second globalization 
strategy proposed in Section 13.41 



3.2 Newton- Algorithm 

The generalized differential dP^a^t] can be defined analogously to (j3.4p and the set-valued 
function dP\^a,b]{ ~ ^PhiUhiv))) contains the characteristic function Xi{v) of the inactive set 

I{v) = {ujen\{- ^PhiVhiv))) (uj) G {a{uj),b{uj)) } . 

By we will denote synonymously the characteristic function Xj(v) ^ well as the self-adjoint 
endomorphism in L^(0) given by the pointwise multiplication with Xx{v)- The Newton-step 
in Algorithm 13.21 now takes the form 

(/ + ^x''SlSh)Sv = -v + P[a,t] (-Iphiyhiv)))^ . (3.5) 

To obtain an impression of the structure of the next iterate = v + 5v we rewrite ()3.5p as 

= P[a,b] (--Ph{yh{v))] - -x^'SiShSv . 

Since the range of 5^ is Y/j, the first addend is continuous and piecewise polynomial (of degree 
k) on a refinement of T/j. The partition is obtained from Th by inserting nodes and 
edges along the boundary between the inactive set T{v) and the according active set, and 
in general contains simplices of higher order than T^. The inserted edges are level sets of 
polynomials of order < k since we assume o, 6 G Y^. 

The second addend, involving the cut-off function x", is also piecewise polynomial of degree 
k on Kh but may jump along the edges not contained in T^. 
Finally f"*" lies in the following finite dimensional subspace of L^(r2) 

K = + (1 - Xl^2 ki, V^2 G Yh] = span {</>,(! - X^U) ■ 

The iterates generated by the Newton-algorithm can be represented exactly with about con- 
stant effort, since the number of inserted nodes varies only mildly from step to step, once 
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the algorithm begins to converge. Furthermore the number of inserted nodes is bounded, see 

MM- 

Since the Newton-increment 6v may have jumps along the borders of both the new and the 
old active and inactive sets, it is advantageous to compute directly, because lies in Y"^. 
To achieve an equation for we add G'^{v)v on both sides of (j3.5p to obtain 

fl+^x''S*r,Sh)v+ = P[a,t] (--Phiyhiv)))) + -x'St.Snv, (3.6) 

and reformulate Algorithm 13.21 as 
Algorithm 3.5 (Newton Algorithm). 

V £ U given. Do until convergence 

Solve (13.61) for v'^ , v := v'^ . 



3.3 Computing the Newton-Step 13.61 

Since defined by (j3.6|) is known on the active set A{v) := \ I{V) it remains to compute 
on the inactive set. So we rewrite (j3.6p in terms of the unknown x^'^^ by splitting u"*" as 

and obtain 

(l + -x''SlSh)x''v+ = Pia,b] ( --Phiyhiv)))) +-x''SlSHV-(l + -x''SlSh){l-xlv+. 
As (1 — x^)^^ is already known, we can restrict the latter equation to the inactive set T{v) 
(x"" + -X^'SlShX^'^v^ = -x'Slz - -x^SlShil - xlv+ . (3.7) 

On the left-hand side of (j3.7p we have now a continuous, selfadjoint Operator on Lp'iX'")^ 
which is positive definite, because it is the restriction of the positive definite Operator 
[l + lx^'SlShX'') to L^l-). 

Hence we are in the position to apply a CG-algorithm to solve (|3.7p . Moreover under the 
assumption of the first iterate lying in 

K\x^ = {x>\v^yh] , 

as does the solution x^v~^ ^ the algorithm does not leave this space because of 



and all CG-iterates lie in Yj^\^^. These considerations lead to the following 
Algorithm 3.6 (Solving ([Ml))- 

Compute the active and inactive sets and . 

Vg G A'' set 

v+{q) = P[„ (^\hiyhiv))iq)) . 
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Solve 



\ a / a a 



,+ 



for x"''J~^ by CG-iteration. By choosing a starting point in Y^^ 
iterates lie inside 17^ L„. 



one ensures that all 



V 



We note that the use of this procedure in Algorithm 13.51 coincides with the active set strategy 
proposed in [3]. 

3.4 Globalization 

Globalization of Algorithm 13.51 may require a damping step of the form 

= V + X{v^ — v) 

with some A > 0. According to the considerations above, we have 



Unless A = 1 the effort of representing vj^ will in general grow with every iteration of the 
algorithm, due to the jumps introduced in each step. This problem can be bypassed by 
focussing on the adjoint state Ph{v) instead of the control v. In fact the function and thus 
also Equation (13. 6p do depend on v only indirectly via the adjoint p = phiy) = S^{ShV — z) 



Now in each iteration the next full-step iterate is computed from (j3.8p . If damping is 
necessary, one computes p'^ = Phiv^) instead of w^. In our (linear) setting the adjoint state 
p~j^ simply is a convex combination of p = Ph{v) and p+ = Ph{v^) 



and unlike the adjoint state p~^ lies in the finite element space Y^. Thus only a set of 
additional nodes according to the jumps of the most recent full-step iterate have to be 
managed, exactly as in the undamped case. 

Algorithm 3.7 (Dampened Newton- Algorithm), v £ U given. 
Do until convergence 

Solve Equation (|3.8p for v~^. 
Compute p+ = ph{yh{v^))- 

Choose the damping-parameter A. (for example by Armijo line search) 

Set p pX = ^P^ + (1 ~ ^)P- 
Algorithm 13.51 is identical to Algorithm 13.71 without damping (A = 1). 

Remark 3.8. The above algorithm is equivalent to a dampened Newton algorithm applied 
to the equation 





(3.8) 



pt = Ap+ + (1 - X)p 




u:=P[aM l-^Ph 
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Another approach, leading a globahzation of Algorithm 13.5^ is to use some globalized, fully 
discrete scheme and then perform 13.51 as a post processing step, compare also [9]. 
Suppose Vh is a discrete approximation to the optimal control u, such that 

\\vh - uWl^u) = C>{h) , 

and let Uhh be its variationally discretized counterpart solving (P^/i). Now, if the q-superlinear 
convergence of order | of the Newton algorithm is mesh independent (see Remark 13. 4p .then 
there exists a radius 6 and a mesh parameter Hq > 0, such that inside the ball Bs{uhh) and 
for h < ho the Newton iteration for Gh converges q-super linearly of order | towards li/^/j. 
Let Uhh be the second iterate of Algorithm 13.51 initialized with Vh. Then, for sufficiently small 
h, we have £ B^{uh) and thus 

9 

\\uhh - u\\l2^p^i < \\uh - Uhh\\L-2{p.) + \\uhh - uWi^^yt) < \\vh - Uhh\\l2(ji) + 0{h^) = 0{h'^) . 

This motivates 

Algorithm 3.9 (Post Processing). 

Solve the fully discretized optimization problem. 
Perform 2 steps of Algorithm 13.61 



3.5 Global Convergence of the undamped Newton Algorithm 

It is not difficult to see, that the fixed-point equation for problem i^hh) 



Uhh = P[ahM (^--^SliShUhh - z)^ 



can be solved by simple fixed-point iteration that converges globally for a > \\Sh\\'j^2(^Q^ L^{n)^ 
see O |6] . A similar global convergence result holds for the undamped Newton algorithm 13.51 

Lemma 3.10. The Newton alaorithm \3.5\ converges globally if a > IHS'p. 

Proof See [13] . I 



4 Numerical examples 



We end this paper by illustrating our theoretical findings by numerical examples. The first 
two examples are solved by Algorithm [321 Algorithm 13. 71 without damping, making use of 
the global convergence property from Lemma [3.10[ The third one involves a small parameter 
a = 10^^ and is hence treated using the globalization strategv 13.71 with Armijo line search. 
Finally the globalization 13.91 is applied at multiple parameters a and mesh parameters h. 
As stopping criterion we require \\P[a,b]{~^P\) ~ ^h\\L^(n) < W~^^ in Algorithm 13.71 using 
the a posteriori bound for admissible v G Uad 

[av + phiy)]^ \i v{uj) = a 
[av+ph{v)]+ if v{uj) = b , 
+ Ph{v) if a < v{lj) < b 



presented in [8] and [10]. 
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Initial Guess 



Step 1 



Step 2 



Step 3 




Example 4.1 (Dirichlet). We consider problem (P) in (II. ip with controls u G L^(0) on the 
unit square Q = (0, 1)^ with a = 0.3 and 6 = 1. Further we set 

z = —Air'^a sin(7r2;) sin(7ry) + (5 o z)r , where r = min (l, max (0.3, 2 sin(7rx) sin(7ry)))) . 

The choice of parameters implies a unique solution u = r to the continuous problem (P). 

Throughout this section, solutions to the state equation are approximated by continuous, 
piecewise linear finite elements on a quasiuniform triangulation with maximal edge length 
h > 0. The meshes are generated through regular refinement starting from the coarsest mesh. 

As discussed in Section [21 problem (P/ih) admits a unique solution Uh and we have 

\\uh - ^^||L2(n) = 0{h^) 
as /i — > 0. There also holds nearly quadratic convergence in L°°(r2) 

\\uh - uWi^iii) = 0{\ log{h)\^h^) 

for domains O C K'^, see [5]. Both convergence rates are observed in Table [H that shows the 
L^- and the L°°-errors together with the corresponding experimental orders of convergence 

^ In ERR{h,^i)- In ERRjhj) 
ln(/ij_i) - ln(/ij) 

for Example 14.11 Lemma 13.101 ensures global convergence of the undamped Algorithm 13.51 
only for a > l/(37r^) ~ 0.0034, but it is still observed for a = 0.001. 

The algorithm is initialized with vq = 0.3. The resulting number of Newton steps as well as 
the value oi C/a for the computed solution are also given in Tabled! 

Figure [1] shows the Newton iterates, active and inactive sets are very well distinguishable, 
the jumps along their frontier can be observed. 

Next we demonstrate another Example, out theory may also be applied to. 



mesh param. h 


ERR 


ERRoo 


EOC 


EOC^ 


Iterations 


Quality 




2.5865e-03 


1.2370e-02 


1.95 


1.79 


4 


2.16e-15 




6.5043e-04 


3.2484e-03 


1.99 


1.93 


4 


2.08e-15 


^/2/64 


1.6090e-04 


8.1167e-04 


2.02 


2.00 


4 


2.03e-15 


^2/128 


4.0844e-05 


2.1056e-04 


1.98 


1.95 


4 


1.99e-15 


^/2/256 


1.0025e-05 


5.3806e-05 


2.03 


1.97 


4 


1.69e-15 


V2/512 


2.5318e-06 


1.3486e-05 


1.99 


2.00 


4 


1.95e-15 



Table 1: L^- and -error development for Example 14.11 (Dirichlet) 
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mesh param. h 


ERR 


ERRoo 


EOC 


EOC^ 


Iterations 


Quality 


V2/I6 


3.9866e-03 


1.1218e-02 


1.94 


1.74 


3 


1.81e-12 


V2/32 


1.0025e-03 


3.2332e-03 


1.99 


1.79 


3 


2.31e-12 


^/2/64 


2.5188e-04 


8.4398e-04 


1.99 


1.94 


3 


9.74e-13 


^/2/128 


6.2936e-05 


2.1856e-04 


2.00 


1.95 


3 


9.37e-13 


^2/256 


1.5740e-05 


5.5223e-05 


2.00 


1.99 


3 


8.91e-13 


\/2/512 


3.9346e-6 


1.3928e-05 


2.00 


2.00 


3 


8.86e-13 



Table 2: Development of the error in Example 14.21 (Nemnann) 



Example 4.2 (Neumann). We next consider an elliptic problem with Neumann boundary 
conditions 

—Ay + y = u in Q , 
dnU = on dVL , 

on 17 = (0, 1)^, with a similar discrete setting as in the previous example. It then is clear, 
how (P) and (Phft) have to be understood. We set a = 1 and choose 

z = — 2(27r^ + l)a cos(7rx) cos(7ry) + (S o i)r , with r = min (l, max ( — 1,2 cos(7rx) cos(7ry))) 
and bounds a = —1 and 6=1. The optimal control to the continuous problem is n = r. 

For a = 1 the undamped iteration still converges globally, although the solution operator has 
norm \\S\\ = 1 as an endomorphism in L'^{Q,). The predicted convergence properties and the 
stopping criterion are the same as above; Algorithm 13.71 is initialized by vq = —1. The first 
four steps of the iteration are displayed in Figure [2] and the behaviour of the approximation 
error between the exact and the semidiscrete solution, as well as the number of iterations and 
the final value of C/c^i is shown in Table El 

The Algorithm has also been implemented successfully for parabolic discontinuous Galerkin 
discretized problems as well as elliptic problems with Lavrentiev-regularized state constraints. 



To demonstrate Algorithm 13.71 with damping we again consider Example 14. H this time with 
a = 10~^. We choose 

MF{p)= p-SlShP[,^,](--p)+S*^z 



as merit function governing the step size of the algorithm. Again we use the same stopping 
criterion as in the previous examples. 



Initial Guess 



Step 1 



Step 2 



Step 3 




Figure 2: The first steps of the Newton-algorithm for Example 14.21 (Neumann) with a = 1. 
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mesh param. h 


ERR 


ERRoo 


EOC 


EOC^ 


Iterations 


V2/2 


1.1230e-01 


3.0654e-01 




_ 


11 


V2/4 


3.8398e-02 


1.4857e-01 


1.55 


1.04 


22 


V2/8 


9.8000e-03 


4.4963e-02 


1.97 


1.72 


19 


V2/I6 


1.7134e-03 


1.2316e-02 


2.52 


1.87 


20 


V2/32 


4.0973e-04 


2.8473e-03 


2.06 


2.11 


33 




8.2719e-05 


6.2580e-04 


2.31 


2.19 


17 


V2/I28 


2.0605e-05 


1.4410e-04 


2.01 


2.12 


20 


V2/256 


4.7280e-06 


4.6075e-05 


2.12 


1.65 


19 


V2/512 


1.1720e-06 


1.0363e-05 


2.01 


2.15 


18 



Table 3: Development of the error in Example 14.11 (Dirichlet) for a = 10 



Table [3] shows errors and the number of iterations for different mesh parameters /i at a 
smoothing parameter a = 10"''. To compare the number of iterations we choose a common 
initial guess uq = 1. The number of iterations appears to be independent of h. 
Finally, to demonstrate the efficiency of Algorithm l3.9t the EOC in the L2(r2)-norm is plotted 
in table HI The disturbances that can be observed for smaller parameter a indicate the decay 
of the environment of q-super linear convergence with decreasing a. 

Acknowledgements 

The first author gratefully acknowledges the support of the DFG Priority Program 1253 
entitled Optimization With Partial Differential Equations. We also thank Andreas Giinther 
for some fruitful discussions. 

References 

[1] K. Deckelnick and M. Hinze. Convergence of a finite element approximation to a state- 
constrained elliptic control problem. SIAM J. Numer. Anal. 45 (5), 1937-1953 (2007). 



h 


1 


10-1 


10-2 


10-3 




10-5 


10-6 


10-7 


10-8 


V2/4 


1.36 


1.37 


1.42 


1.81 


1.81 


1.56 


1.24 


0.96 


-4.07 


V2/8 


1.78 


1.78 


1.75 


1.53 


1.94 


2.41 


2.41 


-6.49 


-2.67 


V2/I6 


1.95 


1.96 


1.97 


2.09 


2.35 


2.68 


2.01 


9.41 


2.96 


y/2/2,2 


2.01 


2.01 


2.00 


1.94 


1.83 


1.78 


2.95 


4.40 


10.43 


V2/64 


1.99 


1.99 


1.99 


2.02 


2.09 


2.19 


2.13 


1.84 


-1.43 


^2/128 


2.00 


2.00 


2.00 


1.99 


1.97 


1.92 


2.01 


2.50 


7.05 


^2/256 


2.00 


2.00 


2.00 


2.01 


2.04 


2.08 


2.11 


2.19 


2.34 



Table 4: EOC of Algorithm 13.91 applied to Example 14.11 (Dirichlet) for different a. 



13 



[2] K. Deckelnick and M. Hinze.: A finite element approximation to elliptic control problems 
in the presence of control and state constraints, Preprint HBAM2007-01, Hamburger 
Beitrage zur Angewandten Mathematik, Universitat Hamburg (2007). 

[3] M. Hintermiiller, K. Ito and K. Kunisch The primal-dual active set method as a semi- 
smooth Newton method. SIAM J. Control and Optim. 13 (3), 865-888 (2003). 

[4] M. Hintermiiller and M. Ulbrich.: A mesh-independence result for semismooth Newton 
methods Mathematical Programming 101, 151-184 (2004). 

[5] M. Hinze: A generalized discretization concept for optimal control problems with control 
constraints, Preprint MATH-NM-02-2003, Institut fr Numerische Mathematik, Technis- 
che Universitat Dresden (2003). 

[6] M. Hinze. A variational discretization concept in control constrained optimization: the 
linear- quadratic case, J. Computational Optimization and Applications 30 (1), 45-61 
(2005). 

[7] M. Hinze, R. Pinnau, M. Ulbrich, and S. Ulbrich. Optimization with PDE constraints. 
Mathematical Modelling: Theory and Applications 23. Dordrecht: Springer (2009). 

[8] K. Krumbiegel and A. Rosch. A new stopping criterion for iterative solvers for control 
constrained optimal control problems. Archives of Control Sciences 18 (1), 17-42 (2008). 

[9] C. Meyer and A. Rosch. Superconvergence properties of optimal control problems. SIAM 
J. Control Optim. 43 (3), 970-985 (2004). 

[10] F. Troltzsch and S. Volkwein: POD a-posteriori error estimates for linear- quadratic 
optimal control problems. Computational Optimization and Applications, Online First 
(2008). 

[11] M. Ulbrich: Semismooth Newton methods for operator equations in function spaces, SIAM 
J. Optim. 13, 805-841 (2003). 

[12] M. Ulbrich:74 new mesh-Independence result for semismooth Newton methods, Oberwol- 
fach Report 4/2009, 78-81 (2009). 

[13] M. Vierling: Ein semiglattes Newtonverfahren fiir semidiskretisierte steuerungs- 
beschrdnkte Optimalsteuerungsprobleme, Diplomarbeit, Universitat Hamburg (2007). 



14 



